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Abstract 

This paper describes and analyzes a method for computing border bases of a zero-dimensional 
ideal /. The criterion used in the computation involves specific commutation polynomials and 
leads to an algorithm and an implementation extending the one provided in This general 

border basis algorithm weakens the monomial ordering requirement for Grobner bases computa- 
tions. It is up to date the most general setting for representing quotient algebras, embedding into 
a single formalism Grobner bases, Macaulay bases and new representation that do not fit into 
the previous categories. With this formalism we show how the syzygies of the border basis are 
generated by commutation relations. We also show that our construction of normal form is stable 
under small perturbations of the ideal, if the number of solutions remains constant. This new 
feature for a symbolic algorithm has a huge impact on the practical efficiency as it is illustrated 
by the experiments on classical benchmark polynomial systems, at the end of the paper. 

Keywords: Multivariate polynomial, quotient algebra, normal form, border basis, root-finding, 
symbolic-numeric computation. 

1 Introduction 

Solving polynomial systems is the cornerstone of many applications in domains such as robotics, 
geometric modeling, signal processing, chromatology, structural molecular biology etc. In these prob- 
lems, the system has, most of the time, finitely many solutions and the equations often appear with 
approximate coefficients. 

From a computational point of view, it is an actual challenge to develop efficient and stable methods 
to solve such problems. First backward stability is expected. The computed solutions should be the 
exact solutions of a system in the neighborhood of the input system. Efficiency is also mandatory 
to tackle the encountered polynomial systems. One can expect for instance that the behavior of the 
method depends mainly on the number of solutions, and partially on other extrinsic parameters such 
as the number of variables. 

To handle the backward stability issue, one may consider classical numerical methods such as 
Newton-like iterations. However these local methods do not provide any guarantee of global conver- 
gence nor a complete description of all the roots. Algebraic methods, on the contrary, handle all the 
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roots simultaneously. They reduce the problem to computing the structure of the quotient algebra 
A of the polynomial ring modulo the ideal / generated by the input system [8]. Such a structure is 
given by a basis B of A as a vector space, and the tables of multiplication in A. Equivalently, it can 
be described by an algorithm of projection of the ring of polynomials K[x] onto the vector space (B) 
generated by B, along the ideal /. We call such a projection, a normal form for the ideal /. 

From the knowledge of the multiplication tables, we deduce either an exact encoding of the roots 
through a rational univariate representation [32] , [12j . or a numerical approximation by eigenvector 
computation, [T], [24j . [34] . Since the eigencomputation can be considered as a numerically well- 
controlled process [13], the challenge becomes now to compute efficiently and in a stable way the 
quotient algebra structure A. 

In the family of algebraic methods, resultant-based techniques (see eg. [T5], [5], [3]) exploit the 
properties of coefficients matrices of monomial multiples of the input equations in some specific degree. 
The table of multiplications are obtained by explicit Schur complements in these matrices [9j . Their 
construction is deeply linked to the geometry of the underlying variety that make these methods very 
tolerant against small perturbations. Unfortunately they heavily rely on genericity hypotheses that 
reduce their applicability. Moreover, the size of the constructed matrices usually grows exponentially 
with the number of variables. 

To avoid such a pitfall, so-called H-bases have also been studied [20] , [22] . They proceed degree by 
degree and stop when the terms of highest degree of the computed polynomials generate the terms of 
highest degree of all the ideal /. The stopping criterion requires the computation of generators of the 
syzygies of these highest degree terms. Though it also yields a basis of the quotient ring A, without 
any apriori knowledge on its dimension, practically speaking, it also suffers from the swelling of the 
size of the linear systems to be solved. 

The approach can be refined further by using a grading for which the highest term of a polynomial 
is a monomial. This leads to Grobner basis computation (see eg. [6]). This approach also yields a 
basis of the quotient algebra A and a normal form for /. As in the other methods, their computation 
can be seen as just a triangulation of a certain matrix. Unfortunately these methods suffer from 
unavoidable instability: the monomial ordering attached to the Grobner basis make the pivot selection 
strategy in the triangulation depend only on the symbolic structure of the rows of the matrix, and 
not on the numerical values of the coefficients appearing in these rows. This can lead to artificial 
unwanted singularities in the representation of the quotient ring (compare for instance the degree 
reverse lexicographic Grobner basis of pi = ax\ + bx\ + e\X\x%, P2 — cx\ + dx\ + E2X\X2 with 
£1 = £2 = and £1 7^ 0, £2 ^ 0). 

To circumvent this artificial difficulty, a new approach based on a normal form criterion which 
involves commutation relations (for bases B connected to 1) was first proposed in [25j. Further 
investigations of this technique [27], [28], including the Ph.D. dissertation [36] lead to a first version [29] 
of a normal form algorithm, which allows to construct efficiently and in a stable way zero-dimensional 
quotient algebra representations. Similar investigations were also pursued in [15j . [16] , but in the more 
restrictive case where the basis is stable by division. Other investigations related to the stabilization 
of the normal form process can also be found in |34] . 

In this paper, we discuss border basis algorithms in the general sense, that is for bases B which 
are connected to 1. See [53] for an introductory presentation of their properties. As in [TT] or [TO], our 
approach is based on linear algebra tools. As for H-bases, we use a grading of the polynomial ring, 
but the construction is optimized in the spirit of [27] . We describes an efficient criterion based on 
commutation polynomials to check the normal form property. This leads to an algorithm, as presented 
in [29], which has been improved to treat optimally the case when a syzygy of the components of 
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highest degree is found, which is not a syzygy of the corresponding polynomials. We prove that the 
syzygies of a (general) border basis are generated by the commutation relations, giving a short and 
concise answer to a conjecture in [15] for basis stable by division. Meanwhile, works related to this 
conjecture for this special case were also investigated in [TJ]. Regarding the numerical stability of 
border bases, we prove that our construction of normal form is stable against small perturbations of 
the input system, if the number of solutions remains constant. 

The paper is structured as follows. In the next section, we recall the notations, used in section [3] 
to prove the stopping criterion for generalized normal forms. In section 01 we show how to recover 
the syzygies from the commutation relations. In section O we recall briefly the maind idea of the 
algorithm described in [2j|] . In section [HI we analyze the stability of this algorithm from a symbolic- 
numeric perspective. Finally, we show the efficiency of our implementation and its numerical behavior 
on classical polynomial benchmarks. 

2 Notations 

We recall some of the definitions stated in [35] , [57] , [35J and add a few more that we will need in the 
sequel. 

Let IK be an effective field. The ring of rt-variate polynomials over K will be denoted by R, 
R = K[x] = K[ari, . . . , x n ]. We consider n-variate polynomials /i, . . . , f s G R- Our goal is to solve the 
system of equations /i = 0, ...,/« = over the algebraic closure K of K. These polynomials generate 
an ideal of K[x] that we call /. The quotient of K[x] modulo I will be denoted by A. From now on, we 
suppose that I is zero dimensional so that A is a finite dimensional IK-vector space. The roots, with 
coordinates in the algebraic closure of K, will be denoted by £i, . . . , Q, with Q = (Q t x, . . . , Ci,n) G K . 

The support supp(p) of a polynomial p G K[x] is the set of monomials appearing in p with non-zero 
coefficients. Given a set S of elements of K[x], we denote by (S) the K- vector space spanned by the 
elements of S. Finally, we denote the set of all the monomials in the variables x = (xi, . . . ,x n ) by 
M. A term is an element of the form A • m with A G K — {0} and m G M. For a subset S of M., 
we will denote by S c the set-theoretical complement of the set S in Ai. For any monomial m G M, 
m ■ M. denotes the set of all monomial multiples of m. 

For any subset S of R, we denote by S + or D(S) the set S + = S U Xi S U ■ • ■ U x n S, dS = S + \S. 

k times 

S + is called the prolongation of S. For any k G N, S>- k ' is S ^ v the result of applying k times the 
operator + on S. By convention, 5^°! is S. 

If B C M. contains 1, for any monomial m & M, there exists k such that m G B^. We say that 
a monomial m is of i?-index k if m G — B^ 1 ^ , and we denote it by 5s(m). 

A set of monomials B is said to be connected to 1, if and only if, for every monomial m in B, 
there exists a finite sequence of variables (^)je[i,q such that 1 G B, VI' G [1, 1], ^j£[i,v\Xij G B and 
n,- e [i,2] x i} = m. 

A set of monomials B is stable by division, if for any m G B and any variable xi such that 
m = Xi m' (m! G M.), we have m! in B. Remark that a set B stable by division, is connected to 1. 

Definition 2.1 Let A be a monoid with a well order relation -<, such that: 

Va, ^,7GA, a^/3^7 + a^7 + /3 

A (A, ^-grading o/K[x] is the decomposition o/K[x] as the direct sum: K[x] = ^ AeA K[x]m, with 
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the following property: 

V/ e K[x] H , g e K[x] [f3] => fge K[x] [a+/3] . 

We denote by degree of f or deg A (/), or deg(/) (when no confusion is possible) and by A(/), the 
following element of A: 

A(/) = min{AeA|/e © K[x] [V] }. 

A'^A 

For any set V C K[x], let V\ = a ,^a K[x] [v] n V. For any A £ A, let A+ = min{A' £ A; K[x]+ C 
K[x] v } and let A" = max{A' £ A; K[x]+ C K[x] A }. 

In order not to be confused with different notions of degree, for any monomial m of K[x] , we define 
the size of m, denoted by m|, to be the integer d such that m = Xi t ■ ■ -Xi d , which itself imposes a 
grading. 

Another classical grading is the one associated with a monomial ordering, where A = N" and -< 
is a monomial order (see [7j p. 328]) such that for all a = (a>i, . . . ,a n ) £ N n , we have K[x][ Q ] = 
Kif 

Definition 2.2 We say that A is a reducing grading if A is a grading and if we have the property: 
for all monomials m,m' £ M, such that m' divides strictly m, A(m') -< A(m), A(m') ^ A(m). 

Both gradings induced by the classical degree and a monomial ordering are reducing grading. Here- 
after, we will denote by A a reducing grading. 

Definition 2.3 A rewriting family F for a monomial set B is a set of polynomials F — {fi}i^x such 
that supp(/i) C B + , fi has exactly one monomial ~f(fi) (also called the leading monomial of f%) in 
9B, ifi(fi) = 7(/j) then i = j, 

Remark that the elements of F can be seen as rewriting rules for the leading monomial using monomials 
of B. 

Definition 2.4 A reducing family F of degree A £ A for a set B is a set of polynomials such that F 
is a rewriting family for B , Vm € dB of degree at most X, 3f £ F \ "/(f) = m. 

For the set B = {1, xq, x\, xqXi}, the set of polynomials F = {.Tq — l,xf — x\,XqXi — x\, xfxQ — x\} 
is a reducing family of degree 3. 

Notice that a reducing family of degree A for a set B (connected to 1) allows to rewrite the 
monomials of (B + )\ modulo F as elements of (B)\. This leads in fact, to the definition of the linear 
projection np, associated to a reducing family for a set B connected to 1. 

Definition 2.5 Given a reducing family F of degree A £ A for a set B connected to I, we define the 
linear projection Ttp : (B + )\ — > (B)\ such that 

Vm £ B\, Trp(m) = m, 

Vm £ dB\, TTp (m) = m — /; 

where f £ F is the unique member of F such as m — "f{f). We extend this construction to (B + )\ by 
¥L-linearity. 
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There is a parallel with the differential algebra terminology, which we want to highlight here. 
To a polynomial f(x\, . . . ,x n ) G R of degree at most S, we can associate the differential equation 
f{di, . . . , d n )$(tx, . . . , t n ) = 0, where di is the derivation with respect to the variable It is a linear 
equation in The solution $(ti, . . . , t n ) lives in the ring K[[ii, . . . , £„]] of formal power series in 
£1, . . . , t n , which we can truncate in degree A >z A(/j). This set J\(t) of truncated series in degree A 
is called the space of X-jets in the literature. See eg. [33J for more details. 

Given a set of polynomials F C K[x] of degree at most A, we consider the so-called solution 
manifold 1Z\ of the differential system F(9)($) = in J\(t) (In our case, it is just a linear space). 
Given A' -< A and $ G J\{t), we can forget the coefficient of degree bigger than A' and project it 
on J\'{t). The set of solutions of F(d)(&) = projected on J\'(t), is then defined by the equations 
F x >(&) = where (F x >) = (F)nI[x] A -. We denote it by w\>(F). We denote by D(F) the new system 
which extends F with the dif(d) = 0, for i G 1 . . . n, f G F, obtained by formal multiplication by the 
di. It is called the prolongation of F. 

The system F of degree at most A is said formally integrable if for any r > 0, tt\(D t+1 (F)) = F. A 
technical condition of involutivity is introduced to ensure the regularity of the differential system [30] . 
A result of Cartan-Kuranishi [U [T71 [3U] asserts that any system F can be transformed by prolongation 
and projection into a (involutive) formally integrable system. The connection of this construction with 
the notion of Mumford regularity of polynomial systems has been detailed in [21] . This correspondence 
has been used explicitly for solving polynomial equations in [31] . This involutive division is also used 
to construct a family of normal form projection related to so-called involutive bases (see eg. [2]). In 
Janet basis construction, however one difference is that the variables do not play a symmetric role. 
So-called multiplicative variables are used to perform reduction and non-multiplicative to extend the 
polynomial vector space. 

In the following, instead of working in the vector space of all polynomials of a given degree, we 
will consider completion procedures based on prolongations and projections relative to a specific set 
of monomials B. Dealing locally with this set of monomials related to the number of solutions of the 
system will allow us to improve significantly the linear algebra stages in this type of methods. 

In the sequel, we will make a heavy use of multiplication operators by one variable that we define 
as follows: 

M itX :(B) x - -> (B)\ 

b I ► TTp(Xib). 

The subscript \ is redundant as soon as we know that F is a reducing family of degree A, and we will 
omit this subscript in the sequel. 

Definition 2.6 Let F = {f\ 1 . . . , f s } be a polynomial set, we denote by F(\^ the vector space: 

F w = ({x a fi\ A(x a fi) < A}). 

Obviously, we have F( X ) C (F)\ where (F) is the ideal generated by F. Next we introduce a definition, 
which is weakening the notion of monomial ordering for Grobner basis: 

Definition 2.7 We say that a function 7 : K[x] — > M (M. is the set of all monomials in the variables 
x\, . . . ,x n ), is a choice function refining the grading A, if for any polynomial p, 7(p) is a monomial 
such that j(p) G supp(p), if m G supp(p), m ^ ~f(p) then j(p) does not divide m, and A(j(p)) = 
max{A(m), m G supp(p)}. The coefficient of the monomial j(p) in p will be denoted by n(p). 
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Example 2.8 In the following, we consider a MacaulajQ choice function 7, such that for all p £ 
K[x], 7(p) = x" 1 • • • x%™ satisfies, \j(p)\ = max{|m|;m G supp(p)} = d, and Bio st- oti — 
maxjcleg^,. (m), m G supp(p) and , \m\ = d\i = 1, ...,n}. In case more t/ian one monomial satisfy 
these conditions, the greatest monomial for the lexicographic order is chosen. 

The monomial returned by the choice function has the same name as the leading monomial of 
an element of a reducing family. This is intended, as we will define a reducing family on the behalf 
of choice functions, and in this framework the two will coincide. Hereafter if S = {pi, . . . ,p s } is a 
polynomial set, then we denote by 7(5*) the set: 7(5*) = {7(^1) ■ • -7(p s )}- 

Definition 2.9 Let 7 be a choice function refining a grading A. For any polynomials p\,P2 G K[x], 
let the C-polynomial relative to 7 and {pi,P2) be 

nl s lcm(7(pi),7(p 2 )) lcm(7(pi),7(p 2 )) 
C(pi,P2) = / w ^ — Pi 7 ^ 1 \ — ^2- 

iei i/ie C-degree of (pi,p 2 ) &e A(lcm(7(pi), 7(p 2 ))) iet i/ie leading monomial of the pair (pi,P2) 
be lcm(7(pi),7(p 2 ))- 

This is almost the same definition as a iS-polynomial [5] when 7 is a monomial ordering. We however 
use a new name to underline that now 7 may not be a monomial ordering (i.e. a total order com- 
patible with monomial multiplication). As we will see in the next section, the C-polynomials express 
commutation conditions for the Mi \ . 



3 Generalized normal form criterion 

Let F = {/1, . . . , f s } be a polynomial system and let I be the ideal generated by F. Remember that 
F(x) (the K- vector space spanned by the monomial multiples of the /j, x a fi of degree ^ A G A) is 
included in I\. Thus, when I\ = Fi\) we can define a normal form modulo /, up to the degree A as 
the projection of K[x]\ along F^ X ) onto a supplementary space (B}\. Hereafter, we consider a set B 
of monomials, containing 1. 

Let F be a rewriting family, and let TL — {m, 3p G F, y(p) = m} be the set of their leading 
monomials then, obviously F allows us to define the projection wf of BUH on B along (F). However 
we may extend this projection using the following extension process: 

Definition 3.1 Let F be a rewriting family. For all m G B, we define ir F (m) = m. For m £ B, 
there exists ml G dB and r integers i\, . . . ,i r G [1, n], such that m = ■ ■ • Xi r m' . We define n F (m) 
by induction on k, as follows. 

• if r = 0, np(m') is defined as np(rn') — 7Tf(w') =m! — f where f G F is such that ~f{f) = m! . 

• Vfe, 1 < k < r, Trp(xi r _ k ■ ■ ■ Xi r m!) — TiF{xi r _ k '^p{xi T _ k+1 ■ ■ ■ Xi r m')), if this latter quantity is defined. 
Otherwise we say that Trp(m) is undefined. 

Remark that the above process allows us to define irp only on monomials, and we extend it implicitly, 
by linearity. Remark also that this extension process is not defined in a unique way. Indeed, two 
different decompositions of a monomial m may lead to two different values of irp (m) . However the 
following theorem shows that this extension process becomes canonical as soon as we check some 
commutativity conditions. 

x It gives monomial basis similar to those given by Macaulay in his multivariate resultant construction |19j . 
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Theorem 3.2 Assume that B is connected to 1. Let F be a rewriting family, and let E be the set 
of monomials m such that for all decomposition of m as a product of variables, 
TrF(xi TTF(xi 1 ■ ■ ■ TTp(xi k ) ■ ■ ■)) is defined. Suppose that for all m £ E and all indexes i , j G [l,rc] such 
that XiXjirn G E , we have: 

TT F (xiTT F (xjm)) = n F (xj7r F (xim)). 

Then tt f coincides with the linear projection its of (E) on (B) along the vector space spanned by the 
polynomials S = {x a f, a G N™, f € F and x a ~f{f) G i?}. 

Proof. Remark that the way we define it, makes tt f inherently a linear multivalued map. Hence to 
prove the theorem we have first to show that under the above hypotheses, np becomes a well defined 
map, and next that this well defined linear map coincide with the projection ns of (E) on (B) along 
(S). 

Remark also that E is obviously stable by monomial division: if all the possible decompositions 
of m as a product of variables m = Xi • • ■ Xi k are such that ir F (xi • • • ir F (xi k )) is defined, then a 
fortiori if m' is a divisor of m, this property is true for m! . Let us show that the extension process 
is independent of the way m is decomposed as a product of variables. Let m = Xi a m' — Xi t m" with 
iq 7^ i\ and m,m',m" G E, then there exists m!" G E (since E is stable by monomial division) such 
that m = Xi x^m'". As m, m! , m", and m'" are in E, Tr F (m'), Tr F (m"), ir F {xi Q m'), irpfe^m"), and 
Tr F (m"') are defined and we have: 

vrf. (x io tt f (m')) = ■Kp(x io Trp(x il -Kp(m'"))), 
ir F (xh n F( m ")) = ^F^ii ^ F (xi Tr F (m"'))). 

The commutation condition guarantees that the two quantities are equal, so that the definition of irp 
does not depend on the way to write m as a product of variables. 

Next we have to show that irp and its coincide on their common set of definition. We do it by 
induction on the size of the monomials: 

It is true that tt|,(1) = tvsO-) = 1 (since 1 G B). For any monomial m ^ 1 in E, the property of 
connectivity of B and the definition of E gives us: 3m' G E and «o G [1, n) such that m = Xi m' and 
np(m') is defined, so that we have: 

7T F (m) - 7r e p(xi m') = def n F (x to n F (m')) =^ucUon ^fa^m')) g (B). 

Now by induction, ml — ns(m') G S\- where A = A(m) and 

m - wp(m) = x io (m' - w s (m')) + (x io irs(m') - np(x ia iis{m'))) G (S). 

Thus np(m) is the projection of m on (B) along (S). □ 
Suppose now that we are given a reducing family of degree 5 instead of a rewriting family. Then 
we can further extend the above theorem with the help of the following lemma. 

Lemma 3.3 Let F be a reducing family of degree A for a set B connected to 1, and suppose that 
V/ G F, A(7(/)) = A(/). With the notation of Theorem ] 3. SI the set E contains the set of monomials 
of degree less than or equal to A . 

Proof. Let m G A4\ be a monomial of degree less or equal to A, then m can be written as m — 
Xi ± . . . Xi d with d — \m\. Let us prove by induction on k < d, that pk = Ttp(xi k Ttp(- ■ ■ 'K F {xi 1 )) ■ ■ ■) is 
defined and that A(pfc) < A(xi k ■ ■ ■ Xi t ). 
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Consider now Xi k+1 pk G B + . For ml G supp(a;j fc+1 pfc) fl B, we have Tr(m') = m' and by induction 
hypothesis A(m') < A(xj fc+1 • • • x^). For m' G supp(x ik+1 pk) fl as F is a reducing family of degree 
A, we have a rewriting rule for to'. The hypothesis that V/ £ F, A( r y(f)) = A(/) implies that to' 
rewrites in terms of monomials of degree bounded by A(a;i fc+1 • • • x^). This proves Pk+i '■= Kp(xi k+1 pk) 
is defined and A(pk+i) < A(xi k+1 ■■■Xi 1 ). 

This proves by induction that irp^x^ ■ ■ ■ TTp{xi d )) is defined, for any decomposition to = x^ ■ ■ ■ Xi d G 
M.\ so that to G F. This ends the proof. □ 

Theorem 3.4 Let F be a reducing family of degree A for a set B connected to 1. If we have: 

• V/gF, A( 7 (/))=A(/). 

• M,, A o M hX - = M a o Mj \- , for l<i,j< n, 

then, we can extend ttf to a linear projection irp from K[x]a onto (B)\ with kernel F(\y 

Proof. As F is a reducing family of degree A, by Lemma I3TB1 we have E D A4\. 

Let us prove that for all to G M.\— and all pairs of indices there exists a way to define up 

such that ir F (xi n F (xjm)) = iTp(xj Trp(xim)). 

As M.\— C M\ C E, np(m) is defined and supp(7r^(m)) C B. We define irp(xim) = 
Trp(xiTrp(m)) and, similarly, Trp(xjm) — 7tF(xj'Kp{in)). With this definition we have: 

■np{x l iTp{x j m)) = M iA (M jtX -(np(m))) 

= M AA (M iiA -(7if.(m))) = irpixjirpixim)). 

which proves the commutation property. We end the proof by applying Theorem 13.21 □ 

Corrolary 3.5 With the hypothesis of Theorem \3.4[ we have K[x]a = (B)\ © -^(A)- 

Let us give here another effective way to check that we have a projection from F/\\ (vector space 
spanned by the monomial multiples of the fi of degree A) onto (B)\ (element of degree A of the vector 
space spanned by B) starting from a reducing family of degree A, without computing explicitly the 
multiplication operators. 

Theorem 3.6 Let A G A. Let F be a reducing family of degree A G A, for B. Assume that V / G 
F, A(7(/)) = A(/) and let irp be the induced projection from (B + )\ onto (B)\. Then V/,/' G F(a) 
such that C(f, f) e(B+) x . 

MC(J,f')) = o 

iff 7Tf extends uniquely as a projection irp from K[x] A onto (B)\ such that ker(-7r|n) = F/\\. 

Proof. By Theorem l3.41 we have to show that this condition is equivalent to the commutation of the 
operators M^a'j A' < A on the monomials of B\-- . 

For any to G B\~- and any i\ ^ Z2 such that Xi 1 m G dB, Xi 2 to G dB, there exists /, /' G F(\-\ such 
that 7(/) = x n to, 7(/') = x l2 to. Thus, we have TT F (x il to) = j(f) - f , TT F (x i2 to) = j(f') - /' and 
C(f, /') = %i 2 f - Xii f S (B)l- Consequently, 

M l2 , x {M luX - (to)) - M iuX (M i2>x - (to)) 

= M i2 ,A(7(/)-/)-M il)A (7(/')-/') 

= 7Tf(^ 2 7(/) - Xi 2 .f) - 7rF(»u7(/0 - x nf) 

= M*iJ' -x l2 f) = n F (C(f'J)). 
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which is zero by hypothesis. A similar proof applies if m £ B or Xi 2 m £ B. 

Conversely, since ker(7rf.) = F^) an d C(f, /') £ F(a) H (-B)^, we have that tt_f(C(/', /)) = 
TTp(C(f, /)) = 0, which proves the equivalence and Theorem 13.61 □ 

Remark 3.7 In this proof, we have shown that if the C -polynomials up to the degree A reduce to 0, 
then the multiplication operators -M^a commute. 

Definition 3.8 A reducing family F for all degrees A £ A on a set B of monomials, connected to 1 
will be called a border basis for B . 

Finally, the previous results leads to a new proof of Theorem 3.1 of [25 : 

Theorem 3.9 Let F be a border basis for a set B of monomials, connected to 1, let ttf be the 
corresponding reduction from (B + ) onto (B) , and let Mi : (B) — > (B) such that V6 £ (B), Mj(6) = 
TTF(xib). Then, 

Mj oMi = Mi o Mj , for 1 < i, j < n 
iff there exists a unique projection irp from K[x] onto (B) such that ker(TT^) = (F) and (^f)\(b+) =^f- 

Proof. Under these hypotheses, by Theorem 13. 4[ for any A £ A, {^f)\(b+) x extends uniquely to a 
projection n F from K[x]a onto (B)\, such that ker(7rf. ) = -F{a)- Since for any A, A' £ A such that 
A ~< A', we have (-Ba')a = B\, and F^ X ) £ {F^x'))\, we also have (np ,)\k[x] x — ^f x - This defines a 
unique linear operator np on K[x] such that 7r[j K [ x ] A = ^p x anc ^ ker(-7r|n) = X^agA ^(A) = (-^)- ^ proves 
the direct implication. The converse implication is immediate. □ 



4 Syzygies and Commutation Relations 

In this section, we analyze more precisely the relations between the polynomials of the border basis 
F = {fu)uedB where f^=uj — p ul with £ (B). These relations or syzygies form a module that we 
denote by 

Syz(F) = {J2 Ke u £ K[x] as ; ^ Kf u = 0}, 

UJ UJ 

where (e u ) ue aB is the canonical basis of K[x] 9B . If F is a border basis family constructed from an 
initial set of polynomials H = {hi, . . . , h s }, we can express / w £ i* 1 in terms of the polynomials in H 
(and conversely) so that a syzygy on F induces a syzygy on H (and conversely). Therefore, we are 
going here to consider only the syzygies on F. 

For any b £ (B) and i £ [l,n], cc, 6 £ can be projected in (S) along F: 

tt f (x 1 b) = x t b - ^2 Mb.uj/uJ- 

More generally, for any p £ K[x], we denote by p l p w the coefficient [i\ w in np(xi b) for b = ttf(p) £ (-B) 
and by p l {p) = YluedB Mp/w Notice that if x,p £ i? then /^ l (p) = 0. 

For any m = Xi 3 ■ ■ ■ Xi k £ _B and ii,i2 £ l...n, the two decompositions ■KF{xi l 'KF{xi 2 m)) = 
■KF{xi 3 -KF{xijn)) yield the syzygy 

x h ^2 ^™ 

)/« = o. (i) 

weas weas weds 
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These relations can also be rewritten as: 



Xi 1 ^' i2 (m) + /Z 1 (xi 2 m) — Xi 2 n 11 (to) — // 2 (x^m) = 0. 

We denote by H the module of K[x] as generated by these syzygies. It is also the module generated 
by the relations of commutation o Mi 2 (b) — Mi 2 o (b) = for all b £ B, ij, 12 £ [1 . . . n]. We 
distinguish the following relations: 

• If Xi 1 m, Xi 2 m £ £?, then y} 1 (m) = fi t2 (m) — 0, [i %1 (xi 2 m) = y}' 2 (xi-\ m) — and to does not 
yield a non-trivial relation of the form |T]). 

• If x^m £ i? but a:i 2 m € dB, and 2;^ to 6 d-B, then 'Kp{xi 2 to) = x^m — f Xi2 m and we have 
the relation 

^ii fxi 2 m ~ fx il Xi 2 m + ^ ] A*x j2 ro,o> ~ 0- (2) 

• If XijTO £ £>, a;i 2 TO £ <9-B, and £i 2 to £ £>, then 7r F{xi 2 m) = x^m — f Xi2 m and we have the 
relation 

2^1 fxi 2 m + ^ ] l J 'Xi 2 m,ufi>> = 0- (3) 

• If Xi x m £ 9_B and x^to £ dB, then ^(x^m) = x^m — fx^mi ^F(xi 2 m) = Xi 2 m — f Xi2 m and 
we have the syzygy 

x i 1 fxi 2 rn ~ Xi 2 f Xii11l — ^ ' (^x^m.^ — ^Xi 2 m,u) f u = ^ (^) 

The syzygies ([5]) and Q are called respectively next-door and across-the-street relations in [16j . The 
syzygies ([3]) does not exist if B is stable under division (since in this case, x l2 m £ B implies X%-y Xj 2 TO £ 
B). All these syzygies are simply the non-trivial relations induced by the "border" C-polynomials 
(see Definition 12. 9[) . 

Let us now prove that the module of syzygies is generated by these commutations syzygies. For 
any monomial m — Xi t ■ ■ ■ x; fc , we can rewrite by induction its projection as: 

7T F (to) = TTp {Xj il 7Tp {Xi 2 •■■)=?«- ^ X n ' ' * ^1-1 ^X^-X^ufw ( 5 ) 

1=1. ..k uiedB 

with the convention that Xi k+1 = 1. We denote by 



^ ] x ii ' ' ' x ii-i ^ ] Mi 
1=1. ..k 



the corresponding element of K[x] as . Notice that this decomposition depends on the order of the 
decomposition m — Xi t ■ ■ ■ Xi k as a product of variables. 

Lemma 4.1 If to = x^ ■ ■ ■ xi k = Xj 1 ■ ■ ■ Xj k then 

^Xi l ---Xi k ^Xj 1 ,...,xj k »— 
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Proof. Consider first a permutation of two variables m — miXi l Xi l+1 rri2 — miXi l+1 Xi l m2 (with 
mi = Xjj • ■ ■Xi l _ 11 ra-i = Xi l+2 ■ ■ -Xi k ). Using ([5]), the two way of projecting m — miXi l Xi l+1 m2 — 
m 1 x ll+1 x ll m 2 yield the syzygy 



mi x 



which is in 5. By iterated permutations of two variables, we transform the sequence Xi x , . . . , Xi k into 
the sequence Xj ly . . . ,Xj k . This allows us to rewrite S Xii ... Xi into ,...,x jk modulo S. □ 

Lemma 4.2 If m9 = m! 9' with 9, 9' G 9-B and m,m' G Ad, then 



2J e w mod S 

uj£dB 



with G K[x] o/ degree < max(|m|, |m'| 



Proof. If m = ■ • ■ Xi d with 9 = Xi d+1 b G <9-B and b £ B, the projection formula ([5]) of m# has the 
form 

TT F (m9)=m9~x ll ---x ld f e ~ x h--- x k-i ^x h+1 -x id e,J^ 



i+i 

1=1. ..d LuedB 

since Xi d+1 6 = and ttf{9) — 9 — fg. If m' = a;^ ■ • ■ , and 0' G 9£? with m9 = m'9', by Lemma I4TT1 
the two decompositions yield the syzygy 

x ii ' ' ' x id e @ ~ x jl ' ' ' x id' e ®' 

+ x ii ' " ' x ii-i 

1=1. ..d w£dB 



x h'" x 3l-i ^x^ ■■■Xj,8',u> e '* 



-l...d' uiGdB 

as an element of H. This syzygy is of the form meg — m' eg' + XLeds e ^ w hh deg(p^) < max(d, d'), 
which proves the lemma. □ 
This yields the following result for a border basis, conjectured in [TH] for the case where B is stable 
by division: 

Theorem 4.3 Let B C M be connected to 1 and let F — (/ w ) w g3b be a border basis for B. Then 
Syz(F) is generated by the relations 0), f3J) and Q). 

Proof. Let a = J2^P^ e ^ e Syz(F). Consider in this sum, any term Xmeg where A 6 K — {0}, 
m e M, 9 G dB and fe(mfl) < |m| + 1. Then, there exists m' G M, 9' G dB such that mO = m'9' 
and 63(171' 9') — \m'\ + 1. This implies that |m| > |m'|. Applying lemma [4~2l the product meg can 
be reduced modulo 2, to 

m'egi + 2J 9^ e w , 
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with deg(q u ) < \m\. Iterating this reduction, we may assume that each term Ame w (A G K — {0}, 
m G supp(p w ), uj G dB) is such that <5b(tow) = \m\ + 1. Notice that the polynomial m/ w has only 
one monomial of maximal £?-index, which is muj. 

As J2luP^ = 0, there exist 9 ^ 9' G <9_B and monomials m G supp(pe), JTi' G supp(pe') such that 
m6 = ml 9' and <5b(to#) is maximal among all terms of the syzygy. By lemma we can replace 
this pair (meg ,m! eg') modulo 3 by a sum of terms of smaller degree. This transformation reduces 
either the number of terms with maximal £?-index or the degree of the polynomials . 

Since we cannot iterate it infinitely, we deduce that a is in 5. □ 

5 Algorithmic issues 

From the preceding sections, we deduce an algorithm as it is done in [29) . The main idea is to translate 
the previous concepts into linear algebra. From Section [3] to compute effectively a normal form, one 
has to find a (monomial) basis B of the quotient algebra connected to 1, and border relations such 
that the multiplication operators commute. The algorithm described in |29j . is a fix point method, 
which updates 

• a potential monomial basis B of the quotient algebra, 

• a set P of polynomials or rewriting rules (with one monomial of their support in dB and the 
remaining monomials in B), 

until a fix point is reached. At each step of the algorithm, the following operations are performed: 

1. The set P' of polynomials of P + with support in B + is computed. 

2. By taking linear combinations, a basis P of the vector space (P'} is computed, such that each 
element of this basis has at most one monomial in dB (and the other in B). 

3. The C-polynomials of the elements of P with their support in B + are computed and reduced 
by P. 

4. If non-zero polynomials with support in B appear, the potential basis B and the polynomial set 
P are updated. The update of B is done by removing some parts of B. The update of P is done 
by combining the elements of P and the reduced C-polynomials, in order to get rewriting rules 
for the new set B. 

The details and the technical proof of termination of the algorithm are available in . We mention 
here that we have improved the algorithm described in this paper since then: the ways the degree 
drops are treated in [29] is now done to avoid to the repeated execution of similar reductions. 

6 Stability of the bases 

This section is devoted to the study of the stability under numerical perturbations, of the bases 
computed by the previous algorithm. 

In many real-life problems, the system f = (/i, . . . ,/ s ) to be solved is given only with limited 
accuracy. However, most of the time, one also knows that the structure of the solutions is invariant 
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in a small neighborhood of the system. Hence one of the feature that is often required to polynomial 
solvers is to produce a representation of the quotient algebra that is stable in a small neighborhood 
of the initial system. The structural numerical stability of the basis is expected in order to have a 
smooth behavior of the coefficients of the representation in the neighborhood of f = (fx, . . . , f s ). By 
definition, a neighborhood of f is an open set in the space of vector of polynomials (hi, . . . , h s ) such 
that A(hj) di A(/j) (i — 1,. . . , n) and which contains f. For e > 0, we define by N e (f) the set of 
systems (hi, . . . , h s ) such that A(hj) ^ A(/,-) and the coefficient vector of hi is at most at distance e 
from the coefficient vector of fi (for the oo-norm). 

Assumption 6.1 Hereafter, 7 denotes a choice function refining a reducing grading A, such that 
for all p G R, j(p) depends only on the support of p and not on the numerical value of its monomial 
coefficients (e.g. Macaulay's choice function, grevlex choice function,. . .). Let^ e be the choice function 
that for any p G R, applies the choice function 7 on the monomials of p, which coefficient norm is 
bigger than e. 

Theorem 6.2 Let f = (/1, . . . , f s ) be a zero dimensional polynomial system such that in a neighbor- 
hood U of f, all systems have the same number D of complex solutions, counted with multiplicities. 
Then for all e > small enough, there exists v > such that for any system f £ N v (f) C U, the basis 
B computed with 7 satisfying Assumvtion \6.1\ for the system f is also a basis for the system f. 

Proof. Let us consider the matrix M whose rows correspond to the coefficient vectors of the 
monomial multiples m fi, with deg(m) + deg(fi) < k + do where k is the number of loops in Algorithm 
[29"] and do the maximum degree of the polynomials /». The columns are indexed by all the monomials 
of degree < n + d . 

We denote by B the monomial set obtained as a basis of A = R/I by applying Algorithm [29] to 
fx,. . . , f s with the choice function 7. By construction, the set of monomials indexing the columns of 
M contains dB. 

Let M 9 be the same matrix as M but constructed with the polynomials fi replaced by generic 
equations, i.e. equations with indeterminate coefficients having the same support as fx, ■ ■ . , f s - Let 
be the block of columns of M indexed by monomials not in B and let N 9 be the corresponding block 
in M 9 . 

Since the operations of Algorithm [29] consist in computing linear combinations of some monomial 
multiples of the polynomials P (see Section [S] and and thus of monomial multiples of fi of degree 
< k + do, the complete computation can be reinterpreted as an optimized triangulation procedure of 
the block N. 

The block N 9 specialized at f is invertible, because any monomial not in B can be reduced by 
the computed border basis of f to an element in (B). This implies that N 9 specialized at f is also 
invertible, for f G A„(f) with v > sufficiently small. Therefore any monomial of dB can be reduced 
modulo the ideal (f) to an element in (B). Consequently, B (which contains 1) is a generating set of 
the quotient algebra A' = R/{?) for any f G N„({). 

As the number of solutions D = \B\ (counted with multiplicity) is left unchanged by small pertur- 
bations in the neighborhood A„(f) of f, the monomial set B, that has exactly cardinality D, is also 
a basis of the quotient algebra A' = R/(f) for the perturbed system f. □ 

Consider now a slight modification of Algorithm [53] : the coefficients whose norm is less than e are 
simply ignored in all the steps of Algorithm [29] . This means that they will not be taken into account 
for choosing a leading monomial, or deciding if a polynomial is nonzero. This small variant will be 
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denoted as the e-algorithm in the next theorem. Remark here that this behavior is quite classical in 
fact, it is more or less what is usually done when neglecting small coefficients in a numerical algorithm. 
Remark also that the following theorem performs rigorously, and that the computed result is not an 
approximation of the true quotient algebra! 
Then the following holds: 

Theorem 6.3 Let f = (/i, . . . , / s ) be a zero dimensional polynomial system such that in a neighbor- 
hood U of f , all systems have the same number D of complex solutions, counted with multiplicities. 
Then for all e > small enough, there exists v > such that for any system f G N u (f) C U, the basis 
B computed with 7 satisfying Assumvtion \6.1\ for the system f is also the basis obtained with 7 e and 
the e-algorithm for the system f. 

Proof. By Theorem 16.21 B is also a basis of the quotient algebra R/(f). This basis B is obtained 
by applying Algorithm 29J to f . The result of this algorithm does not change for the system f if we 
replace the choice function 7 by 7 € for e > small enough (eg. smaller than the minimum of the norm 
of the coefficients of the polynomials on which 7 is applied) . 

Let us show by induction on the loop index k of the algorithm that the steps and the polynomials 
computed by the e-algorithm with 7 e are the same as for the direct algorithm, up to the terms of norm 
smaller that e. 

It is true obviously for the first step k = 1. Let us suppose now that steps 1, . . . , k' of Algorithm 
[29] ran on f and its e variant ran on f ' are structurally the same and let us show that step k' + 1 is also 
structurally the same for the two computations. The coefficients of all these constructed polynomials 
are rational functions of the coefficients of f, which are well defined in a neighborhood of f . If f is 
close enough to f , the monomials which are in the support of the constructed polynomials for f but 
not in the support of the constructed polynomials for f , have coefficients of norm smaller then e. 

If, by hypothesis, the first k' steps are structurally identical, the same polynomials, up to terms of 
norm smaller than e, appear when selecting the polynomials in P + , and the same C-polynomials are 
constructed (see Section [5] and [H]). By continuity of the coefficients of the constructed polynomials, 
the same pivots (of norm bigger than e) are used to construct the new elements in P. Similarly, in a 
neighborhood of f , up to terms of small norm, the C-polynomials not reducing to zero are the same for 
the two computations. If choices of leading monomials are to be performed, then, by Definition l6.ll "f e 
will select the same monomials for f ' and f . Finally at the end of step k' + 1 the same computations 
are performed, up to terms of small norm and the coefficients of the new constructed polynomials are 
rational functions of the coefficients of f, which are well defined in a neighborhood of f. 

This ends the induction showing that the two computations are structurally identical for f G AT„(f) 
with v > small enough. Hence B will also be found as a basis of R/(C) with this e-algorithm. □ 

The rewriting rules obtained from the e-algorithm are close to the exact rewriting rules of the 
system f. Their numerical quality can be improved by iterative refinements such as Newton-like 
iterations using the commutation relations. Such approach has been investigated in [34) . 

Remark 6.4 (Numerical certification) Theorem 16.51 shows the continuity of the normal form 
computation with respect to the coefficient of the input system. It states that there exists a region 
of stability for the computed quotient algebra representation, but it is an open problem to compute 
apriori the value e and v for a given polynomial system in order to control the size of the allowed 
perturbations. This problem is the subject of further work. 

Remark 6.5 (Flatness) Theorem \6. 3\ also shows that if we consider a rationally parametrized family 
of systems ft G N u (f) for all t G [0, 1], such that fo = f and Fq is the border basis for B, then the 
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set B is also a basis of At = R/(it)- Moreover, the border basis F t of ft for B is of the form 
F t = (lu — Pu,t)u£dB> where p u ,t S (B) is a continuous (rational) function oft on [0,1] such that 
Fo = (to — po,u>)u>edB- This also implies that the C -relations (0j and ^ generating the syzygies of F t 
are continuous (rational) functions oftE [0, 1] , which coincide with the C -polynomials relations of Fq 
at t = 0. Consequently any syzygy of Fq which is a combination of the C -relations can be deformed 
continuously into a syzygy of F t (for t £ [0, 1]/ In other words, the systems f in the neighborhood 
N v (f) are flat deformations of the system f [7|/. 

7 Experimentation 

The algorithm described in the previous section is implemented in the library SYNAPfQ. It corre- 
sponds to about 50 000 lines of C++-code. It involves a direct sparse matrix solver. The numerical 
approximation of the roots are obtained by eigenvalues computation, using the library lapack (the 
routine zgees) and the strategy described in [5j. The computations are performed on an AMD- Athlon 
2400+ with 256MB of main memory. We show the results obtained with our implementation in the 
case where the grading that we use for K[x] is the usual one. In the sequel, drvl will refer to the choice 
function associated to the Degree Reverse Lexicographical order, dlex to the degree lexicographical 
order, Mac to Macaulay's choice function (see Example 12 . 8[> . minsz to the choice function over the 
rational that minimizes the memory needed in the reduction loop (this choice only minimizes a local 
step and does not insure local minimality of the global required memory), and mix to the choice 
function that returns randomly either the result of minsz, or of drvl applied to its input. To anal- 
yse the quality of approximation, we mesure the maximal norm at the computed roots of the initial 
polynomials /j and denote it hereafter by mnacr. 

7.1 Generic equations 

The method we propose here is an extension of the Grobner bases computations. As such it can 
compute Grobner bases. The implementation we have is not as optimized as the Grobner bases ones 
that are being worked on for decades. An important work, mostly on linear algebra, remains to be 
done on our program. However we want to show that the method we propose here is competitive, 
and that it does not lose the good practical efficiency of Grobner bases computations [TU]. As the 
arithmetic used in our programs for doing exact computation is the rational arithmetic of GMP, which 
is much slower than integer computations used in the other software we will restrict our-self to the use 
of modular arithmetic. The family of examples we have chosen is the Katsura(n jf| equations. These 
equations are projective complete intersection with no zero at infinity. Using the Macaulay choice 
function, we know apriori that Macaulay's basis will be a monomial basis of the quotient algebra, so 
we know, apriori, what monomials will be leading monomials for the whole computation; so in this 
case we can guarantee that no test to returns erroneous result even using floating point arithmetic. 
We compare first our program to one of the best implementations available, Magma's implementation 
of F4 algorithm [TU] . 

2 http: //www-sop. inria. fr/galaad/synaps/ 
^http: //www-sop. inria. fr/galaad/data/ 
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n 


Synaps mac 


Synaps drvl 


Magma drvl 


7 


0.19s 


3M 


0.22s 


3M 


0.05s 


3M 


9 


6.17s 


5M 


8.44s 


5M 


1.670s 


7M 


10 


32.39s 


14M 


56.84s 


13M 


13.50s 


23M 


11 


252.05s 


50M 


387.97s 


45M 


96.76s 


70M 


12 


1935.25s 


191M 


3072.08s 


157M 


1560.76s 


240M 



Let us mention that Gb, one of the reference implementation of Buchberger's algorithm, spends 
659s on Katsura{ 10). 

Numerically we observe that choosing the mac function also results in a better conditioning of the 
computations. More precisely on Katsura(6), and using a threshold of 10 -10 we have: 



7 


drvl 


dlex 


mac 


drvl 


dlex 


mac 


drvl 


dlex 


mac 


# bits 


128 


128 


128 


80 


80 


80 


64 


64 


64 


time 


1.98s 


2.62s 


1.64s 


1.35s 


3.98s 


0.95s 






0.9s 


mnar 


1Q -^» 


w-' M 


10 -au 


10 -au 


10" ib 


10 -iu 






io- n 



For the 64 bits computation the results computed for the drvl and dlex orders are erroneous due to 
roundoff errors. The time given is the time spent in the computation of the multiplication matrices. 
Afterward, we used either LAPACK to perform the eigenvector computations or Maple when we needed 
extended precision. Because of the different nature of these tools, we do not report on the solving 
part timing. Finally we show here the amount of memory needed to perform the computations over 
Q, using GMP mpq. 





mac 


minsz 


drvl 


mix 


time 


4.22s 


30.21s 


6.54s 


7.83s 


size 


4.2M 


6.1M 


4.4M 


4.9M 



On these experiments, we observe that the local strategy minsz which tends to minimize locally 
the size of the coefficients in the linear algebra operations, does not yield globally the optimal output 
size. In this example, the time and the memory size seem to be correlated. 

7.2 Parallel robot 

Let us consider the famous direct kinematic problem of the parallel robo10 [53] . First we use floating 
point numbers to check the numerical requirements of the computations for different orders. For 
testing a number to be 0, we will use a leveling (here 10 -8 is enough) and we will check afterward 
that the choices performed are the same as those done using modular arithmetic. This is equivalent 
to the use of an hybrid arithmetic [35] . 



7 


# bits 


time 


mnacr 


drvl 


128 


2.07s 


0.3 * 10~ 24 


dlex 


128 


4.27s 


0.3 * 10~ 2;i 


mac 


128 


2.22s 


0.1 * 10 -24 



Here we see that choosing the right choice function can increase (but not so much in this case) the 
numerical accuracy of the roots. Hereafter we use the parametrization of [18] for solving, it involves 
more variables, gives better timings but less correct digits on the final result. 

4 http: //www-sop. inria. fr/galaad/data/ 
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# bits 


time 


mnacr 


250 


1.32s 


1Q -B3 


500 


2.23s 


10 -14(J 



Finally we performed tests using rational arithmetic. 



7 


mac 


minsz 


drvl 


mix 


time 


315s 


229.08s 


201.65s 


257.50 


size 


17M 


14M 


16M 


13M 



In fact, it is not so surprising to see that the choice function 7 has a big impact in terms of the 
computational time and of the memory required. However in this problem, the time and the memory 
size do not seem to be correlated as in the previous case. 

We also mention here that over-constraining the system can result in a dramatic decrease of the 
computation time. Indeed expressing more constraints than necessary can simplify computations 
significantly (see [29]). 

Acknowledgments: We thank A. Quadrat for interesting discussions on differential algebra, prolon- 
gation and involutivity. 
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